In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import curve_fit
Models
In [100]:
def gauss(x, a, center, sigma, offset):
return a * np.exp(-(x-center)**2/sigma**2) + offset
def twogauss(x, a1, center1, sigma1, offset1, a2, center2, sigma2, offset2):
return gauss(x, a1, center1, sigma1, offset1) + gauss(x, a2, center2, sigma2, offset2)
def gausslinear(x, a, center, sigma, offset, slope, intercept):
return gauss(x, a, center, sigma, offset) + slope * x + intercept
def twogausslinear(x, a, center, sigma, offset, a1, center1, sigma1, offset1, slope, intercept):
return gausslinear(x, a, center, sigma, offset, slope, intercept) + gauss(x, a1, center1, sigma1, offset1)
In [34]:
number_of_events = 1000
bin_range = [-10, 10]
bin_width = 0.3
In [35]:
bins = np.arange(bin_range[0], bin_range[1], bin_width)
bin_mids = np.arange(bin_range[0], bin_range[1], bin_width) + bin_width/2.
In [69]:
line_strength = [0.4, 0.6]
line = np.random.randn(number_of_events*line_strength[0])
background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[1])
photons = np.append(line, background)
In [70]:
hist, bin_edges = np.histogram(photons, bins = bins)
In [71]:
ydata = hist
xdata = bin_mids[:-1]
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.plot(bin_mids[:-1], hist)
Out[71]:
In [72]:
ydata = hist
xdata = bin_mids[:-1]
fit_guess = [1, 0, 1, 0, 1, 1]
popt, pcov = curve_fit(gausslinear, xdata, ydata, sigma=np.sqrt(ydata)+1, p0=fit_guess)
In [73]:
popt
Out[73]:
In [77]:
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.ylim(0)
#plt.plot(bin_mids[:-1], hist)
plt.plot(xdata, gausslinear(xdata, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]), label='Fit')
plt.legend()
Out[77]:
In [75]:
perr = np.sqrt(np.diag(pcov))
In [76]:
perr
Out[76]:
In [101]:
number_of_trials = 10000
centers = np.zeros(number_of_trials)
centers_errors = np.zeros(number_of_trials)
def random_gauss_fits(number_of_trials):
line_strength = [0.4, 0.6]
for i in np.arange(0, number_of_trials):
line = np.random.randn(number_of_events*line_strength[0])
background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[1])
photons = np.append(line, background)
hist, bin_edges = np.histogram(photons, bins = bins)
fit_guess = [1, 0, 1, 0, 1, 1]
popt, pcov = curve_fit(gausslinear, bin_mids[:-1], hist, sigma=np.sqrt(ydata)+1, p0=fit_guess)
centers[i] = popt[1]
return centers
In [102]:
centers = random_gauss_fits(number_of_trials)
In [110]:
plt.hist(centers, bins=20, label='$\sigma$ = ' + "%0.3f" % centers.std(), range=[-0.5,0.5])
plt.xlabel("Fit center")
plt.xlim(-0.5, 0.5)
plt.legend()
Out[110]:
In [130]:
number_of_events = 1000
line_strength = np.array([0.4, 0.4, 0.2])
line_positions = [0, 10]
bin_range = [-10, 20]
bin_width = 0.3
In [131]:
bins = np.arange(bin_range[0], bin_range[1], bin_width)
bin_mids = np.arange(bin_range[0], bin_range[1], bin_width) + bin_width/2.
In [132]:
line1 = np.random.randn(number_of_events*line_strength[0]) + line_positions[0]
line2 = np.random.randn(number_of_events*line_strength[1]) + line_positions[1]
background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[2])
photons = np.append(line1, np.append(line2, background))
In [133]:
hist, bin_edges = np.histogram(photons, bins = bins)
In [134]:
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.plot(bin_mids[:-1], hist)
Out[134]:
In [10]:
popt, pcov = curve_fit(twogausslinear, xdata, ydata, sigma=np.sqrt(ydata)+1, p0=[1, 0, 1, 0, 1, 10, 1, 0])
In [11]:
popt
Out[11]:
In [12]:
pcov
Out[12]:
In [13]:
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.ylim(0)
plt.plot(bin_mids[:-1], hist)
plt.plot(xdata, func(xdata, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], popt[7]), label='Fit')
plt.legend()
Out[13]:
In [17]:
pcov
Out[17]:
In [14]:
perr = np.sqrt(np.diag(pcov))
In [15]:
perr
Out[15]:
In [16]:
x =
plt.plot()
In [ ]: